Interface dynamics of a two-component Bose-Einstein condensate driven by an external force 
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The dynamics of an interface in a two-component Bose-Einstein condensate driven by a spatially uniform 
time-dependent force is studied. Starting from the Gross-Pitaevskii Lagrangian, the dispersion relation for 
linear waves and instabilities at the interface is derived by means of a variational approach. A number of diverse 
dynamical effects for different types of the driving force is demonstrated, which includes the Rayleigh-Taylor 
instability for a constant force, the Richtmyer-Meshkov instability for a pulse force, dynamic stabilization of the 
Rayleigh-Taylor instability and onset of the parametric instability for an oscillating force. Gaussian Markovian 
and non-Markovian stochastic forces are also considered. It is found that the Markovian stochastic force does not 
produce any average effect on the dynamics of the interface, while the non-Markovian force leads to exponential 
perturbation growth. 
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I. INTRODUCTION 



Superfluid dynamics of interfaces between Bose-Einstein 
condensates (BECs) of dilute atomic gases has gained active 
interest, both in experimental lUl-llt], and theoretical works |0- 
[13I1 . These studies started with the experimental realization 
of a multicomponent BEC (Ulltl and the theoretical investiga- 
tion of the planar interface structure ||^-^ . At present, a good 
deal of the research focuses on the hydrodynamic instabili- 
ties at the interfaces in BECs, such as the Kelvin-Helmholtz 
instability, the Rayleigh-Taylor (RT) instability, the ferrofluid 
normal-field instability, and the Richtmeyer-Meshkov (RM) 
instability lfl2l - [T7ll . A common feature of these fundamen- 
tal instabilities in the classical gas dynamics is the produc- 
tion of vortices at the unstable interface. In that sense, dy- 
namics of an interface between two quantum fluids has some 
unique features with quantized vorticity being, probably, the 
most specific and the most interesting among them. At the 
same time, there is also much similarity between the classical 
and quantum gas dynamics, which allows "borrowing" some 
classical results and applying them directly or with modifi- 
cations to the quantum case, see Refs. il2l4l5ll . Such "bor- 
rowing" is especially typical in studies of the linear instabil- 
ity stages, for which the difference between the classical and 
quantum cases is expected to be minimal. Still, the problem 
of rigorous derivation of such effects in BEC systems from 
the basic principles of quantum theory remains open for most 
cases. The derivation is needed not only to justify the ex- 
trapolations from the classical to quantum gas dynamics, but 
also to find limits of such extrapolation as well as to indi- 
cate intrinsic quantum effects involved into the problem. As 
an example, we point out the study of the interface waves in 
a system of two BECs undertaken in Ref. |8]. Apart from 
common capillary waves anticipated from the classical gas- 
dynamics and related to in-phase perturbations of two BECs, 
Ref. 1 8] demonstrated also the possibility of quantum counter- 
phase perturbations, for which the local interpenetration depth 
of two wave functions for different BECs oscillates in time. 
Thus, one purpose of the present paper is to provide rigor- 
ous derivation of the dispersion relations for some fundamen- 
tal quasi-hydrodynamic instabilities in a two-component BEC 
driven by a time-dependent force. 



As a particular experimental realization of such a force, 
Sasaki et al. 1 12] suggested a system of two phase-segregated, 
interacting BECs with different spins placed in an external 
magnetic field gradient. This is possible, for example, for 
two magnetic sublevels of the F — \ hyperfine state of ^^Rb 
with \F,mf) ~ 1 1,-1) and Starting from the Gross- 

Pitaevskii (GP) action for the two component BEC system, 
we derive the equations of motion of small perturbations of 
the interface. 

The other purpose of the present paper is to demonstrate the 
diversity of hydrodynamic effects in a system of two separated 
BECs under the action of a time-dependent force. Recently, 
within classical gas dynamics, there has been increased inter- 
est in the problems of the modified RT and RM instabilities 
subject to time-dependent gravity, e.g. see lfT8l - l2Tll . The con- 
figuration of time-dependent gravity was discussed in the con- 
text of specially designed laboratory experiments [22.1. in the 
scope of the problem of inertial confined fusion lll9ll23L l24ll 
and in relation to flame dynamics ll25l - l28ll . In particular, stabi- 
lization of the RT instability by an oscillating high-frequency 
addition to the gravity acceleration has been proposed both 
in the traditional RT configuration |29] and for the inertial 
confined fusion 1231 l24ll . Development of the parametric in- 
stability at a flame front under the action of acoustic waves 
and possible stabilization of the hydrodynamic flame instabil- 
ity has been investigated experimentally and theoretically in 
Refs. ll30l - t32ll . Here we show that similar effects take place 
in BECs. Apart from the RT and RM instabilities considered 
before in llT2l[l5ll . we demonstrate also the possibility of dy- 
namical stabilization of the RT instability by high-frequency 
oscillations, as well as triggering of the parametric instability. 
Another interesting option, which has not been considered yet 
even in the classical gas-dynamics corresponds to a Gaussian 
stochastic force, which may be Markovian or non-Markovian 
with zero or non-zero time correlations, respectively. In the 
present paper we show that the Markovian stochastic force 
does not produce any average effect on dynamics of the in- 
terface, while the non-Markovian force leads to exponential 
growth of perturbations. 

The paper is organized as follows. In section II we for- 
mulate the model. In section III we consider a planar (one- 
dimensional) quantum system with a magnetic gradient push- 
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ing the BECs towards each other, for the case of both constant 
and time-dependent external forces. In section IV we derive 
the equation, which describes the dynamics of small interface 
perturbations. In section V we consider interface dynamics 
under the action of different forces, including a constant force, 
a harmonic force, a pulse force, and a noisy Gaussian force. 
The results of the paper are summarized in section VI. 



II. THE MODEL 

The system of two BECs at zero temperature is described 
by two macroscopic wave functions ^j(f, r), and the mean- 
field GP Lagrangian [33]: 
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2 y j dt 

-y;(f,r)|*/-^|*,-|' 




(1) 



where Vj{i, f) are external potential energies of the BEC com- 
ponents. The tildes in Eq. ^ designate dimensional vari- 
ables, which will be scaled in the following. Keeping in mind 
a system with two magnetic sublevels of ^^Rb, we assume the 
same particle masses for both components mi =m2 = m and 
the same interaction parameters i = g22 = g- The interaction 
parameters gn, g22, gn are related to the respective s-wave 
scattering lengths atj by gij = A-TtfP'ajj/m. A dimensionless 
interaction parameter is introduced through the definition 



r=gi2/g-i- 



(2) 



The two components are separated in space if §12 > ^/giigii, 
which implies 7 > 0, see |33]. Throughout this paper we use 
the condition of weak separation 7<C 1. Modification of all 
results of the present paper to the case of different BECs is 
straightforward. 

The unperturbed flat interface between two components of 
BEC is located in the {x,y) plane with z — indicating the 
plane of symmetry. Fig. 1 . The components 1 and 2 are char- 
acterized by their projection of the atom magnetic moment on 
the jc-axis nip = ±1, respectively. 

Using the states |F,m/r) 1) and |1,1) of the atoms 

makes it possible to apply the Stern-Gerlach force, push- 
ing the components in the opposite directions perpendic- 
ularly to the magnetic moment quantization axis. We 
use a magnetic field whose magnitude has the gradient 
V|B(f)| = xdBx{t)/dz = i(.B'{t), which is uniform in space but 
may be time-dependent. For the case of ^^Rb atoms the Stern- 
Gerlach potential energy in Eq.© reads 



v,(f..)=j^^<<2(-.)', 



(3) 



where fls is the Bohr magneton, and 1/2 is the Lande factor. 

The chemical potentials fij in Eq.([T]i define normalization 
of stationary wave functions ^'j, and are equal to energy per 




FIG. 1: Schematic of the two-component BEC in an external non- 
uniform magnetic field. 



particle of j-th condensate, which is made up of the interac- 
tion energy due to s-wave scattering and of energy of the par- 
ticle in the external potential. In order to study low-energetic 
hydrodynamic processes in BECs we assume that energy per 
particle is changed negligible due to the external potential, and 
for normalization one can use value of the chemical potential 
for B' = in an infinite system. Such system is characterized 
by a uniform density sufficiently far from the interface (i.e. on 
length scales much larger than the interface thickness); then 
the chemical potential is given by a usual result for a uniform 
BEC with number density no, jJ. ~ gfiQ. However, real BECs 
are always finite; if the length scale of our system along z-axis 
is 2Lj (each condensate is of size L,), then the condition that 
the external potential is a small perturbation, reads 



Lz — < gno- 



(4) 



Dimensionless units are introduced by scaling coordinates 
by h/y^mjj., time by Ti/jj., and atom number density (or con- 
centration) by n/g. The dimensionless magnetic force acting 
on the components is 



bi (t) 



~b2{t) 



TijXBB' (f ) 



-b{t). 



The dimensionless Gross-Pitaevskii Lagrangian for two 
macroscopic wave functions {t , r) at zero temperature reads 



E 

.7=1.2 



Vi'A 



2 I ^ dt 



<5^* 
11/ . J 

' dt 



-(l+7)|>I'i|'|>I'2r 



(5) 
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and equations of motion are obtained from the action principle for z> Zp. With these approximations Eq.(l9|l becomes: 



SS _ d 

5^* ~ dt 



SS 



{^^'*/^t) 



(6) 



where the action is S = J dtJ^{t). The result is the coupled 
GP equations: 



---i-z/.(o+i^ir+(i+7)i'P2r 

-^-l+zK0 + l^2|' + (l + 7)l^i|' 



'Pi, (7) 
^2. (8) 



III. 



PLANAR DENSITY PROFILES FOR A 
TWO-COMPONENT BEC 



The present work is devoted to multidimensional instabili- 
ties bending an interface between two BEC components in an 
external magnetic field gradient. As the first step in the study, 
we have to specify the planar (one-dimensional, ID) quasista- 
tionary state, which is influenced by a spatially uniform exter- 
nal field. In the subsections A and B we consider the cases of 
constant and time-dependent external fields, respectively. 



A. Planar profiles in the case of a stationary force 

Here we study the ID hydrostatic state of BECs, when the 
constant external force bo is balanced by internal forces of the 
condensates. Let us assume that such a state was created by 
adiabatic switching on of the external field, and no perturba- 
tions bend the interface. The steady state is described by two 
real- valued functions 1/0^(2), and the system of dimensionless 
GP equations reads 

0== + ™-(l + 7)ro2roi, (9) 

p ^ 1 t/V|2 ^ _ ^Ij^^^^ (\j^y) ^Q2. (10) 

2 az 

The time-independent force pushes the BEC components to- 
wards each other when bQ > 0. With our choice of equal 
masses and scattering lengths for both components, the prob- 
lem possesses the symmetry i//oi (z) — ij/oi i^z). In the case 
of zero magnetic gradient the density profiles have been found 
in 1^ within the hmit of weak separation 7 <C 1 as 



noi (z) = "02 (-z) 



1 



exp 



-1 



(11) 



We calculate the wavefunction of condensate 1 in the bulk 
of condensate 2. According to Eq. (|9]l, the wave functions 
penetrate into each other with the characteristic depth Zp = 
l/^/Ty, see also (tI], which means that Xj/oi ^ 1 for ^ > ^p- 
The density profile of the bulk is described by the Thomas- 
Fermi approximation (TEA), which means that Xj/m ~ 1 — 



1 Vol 

2 dz^ 



[7-(2 + 7)/7oz]roi- 



(12) 



The type of solution to Eq.(fT2b changes from exponentially 
decaying to oscillating in space, as z becomes larger than 
Zo = 7/ [(2 + Y)bQ] w 7/(2/70)- However, the magnitude of the 
spatial oscillations depends on the value y/Qi (Zo)'- it decreases 
exponentially with increasing Z^. Therefore, the spatial os- 
cillations become signific ant w hen Zo becomes comparable 



with Zp, or when = 7\ 



in this case the TEA becomes 



invalid, and Eqs.(l9]l.(fT0ll must be solved exactly. The spa- 
tial oscillations may be interpreted as quantum interference of 
matter waves, caused by interplay between the non-linearity 
and the dispersion (the kinetic term). This phenomenon is be- 
yond the scope of the present work, and it will be studied in 
detail elsewhere. Thus, here we use parameters satisfying the 
condition 



bo < bocr = rVr/^, 

or the magnetic gradient B' satisfying 



(13) 



(14) 



which means that the spatial oscillations of the hydrostatic 
profiles may be neglected. It follows from Eq. ( fT4l l that 
the critical magnetic gradient Bq^.,. depends on the differ- 
ence between inter- and intra-species interaction parameters 
7 = {gn^ S)ls^ on density via the chemical potential 
since /i « guQ in the TEA for an infinite system in a suffi- 
ciently weak external field. Therefore, the effect of quantum 
interference becomes stronger at lower densities of the BECs. 
Equation ( fTST i yields the dimensionless condition of a well- 
defined interface between two condensates subject to the par- 
ticular external force. In the system of two magnetic sublevels 
of ^^Rb we have 7~ 10 ^ (see 12]), and the typical density is 
5 • 10''*cm^^'. The numerical solution to Eqs. (|9]l, (fTOl i for 
these parameters shows that the spatial oscillations are notice- 
able at bo ~ lO^"*; so we take b < 10^'*. Numerical solutions 
for the wave functions near the interface are presented in Eig. 
2 for bo = 7. 5- 10^^; 10^^, plots (a) and (b), respectively. The 
dashed lines in Eig. 1(a) show the approximate analytical so- 
lution of Eq.lfTTI) obtained for the case of zero magnetic field. 
As we can see, the analytical solution Eq.lfTTTi fits the numer- 
ical results quite well even in the case of non-zero magnetic 
field. 



B. ID Bogoliubov excitations of BEC due to a harmonic force 

Since we are interested in the development of hydrody- 
namic instabilities for time-dependent external forces, we 
need to know how a harmonic force influences the planar 
density profiles of BECs. Similar to the constant force stud- 
ied above, for time-dependent forces we also consider rather 
small magnetic fields, so that the respective ID Bogoliubov 
excitations may be interpreted as linear acoustic waves. An- 
other type of a time-dependent force studied in this paper is a 
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FIG. 2: Hydrostatic profiles of two components pushed to each other 
by a constant force, b = 7.5-10^^; 10^^, figures (a) and (b), respec- 
tively. The dashed lines show the analytical solution Eq. ( I14t . The 
inset in figure (b) presents the Airy function of the second kind. 



pulse force corresponding to a weak shock in hydrodynamics. 
Still, it is known from hydrodynamics that a weak shock may 
be also described as a hnear acoustic wave llsill . For this rea- 
son, in the present subsection we study modifications of the 
BEC density in the ID geometry due to the harmonic force. 

Because of the problem symmetry with respect to z = 0, we 
consider a one-component semi-infinite BEC for z < con- 
fined by a wall at z = 0. The magnitude and frequency of the 
force determines the energy of the generated Bogoliubov exci- 
tations. In this subsection we express density and velocity per- 
turbations, 8n and 5U, in terms of magnitude and frequency 
of the external force, and find the restrictions on the driving 
force, for which the density perturbations are small 5n <^ no- 

In the case of a harmonic force, the dimensionless external 
potential is z^.sexp(/£lf), where bs is the force amplitude and 
Q. is the frequency. Our model system is then described by a 
time-dependent GP equation 



.(91// 1 d^\(/ 



dt 



2 dz^ 



W*W^ + WVtrap (z) +zb,exp{iD.t)\ir. (15) 



where V,rap (z) 



+00 for z < and V,rap (z) = for z > 0. 
Since the trap potential Vnap (z) experiences discontinuity at 



the "wall", then the wave function and density are specified 
only for z > 0, no healing of the condensate wave function is 
required near the wall in TEA, and equilibrium is specified by 
the uniform density no ~ I- In that case the velocity pertur- 
bation 5U{z,t) is zero at the wall 5U (0,f ) — 0. The perturba- 
tions of density and velocity, dn and dU, may be expressed in 
terms of the perturbation wave function 5y 



V/(z,f) =exp(-/f)(Vno + 5V^)> 



(16) 



where 5\j/{z,t) — u (z) exp {iHt) + v* (z) exp(— /fif), and m(z), 
V* (z) are time-independent functions. Since the perturbations 
are small, we find 

5n — exp(inf) (m + v) +C.C., 



OU ~ -zr- — — exp(;£2f) — {u — v) —c.c, 
oz 2i dz 

where c.c. denotes complex conjugate. Then the boundary 
condition at the wall reads 



dz 



(m - v) 



= 0. 



(17) 



We linearize the GP equation ( fTSl ) 

.d ^ 1 ^2 

and obtain the Bogoliubov equations 



8\jf + 25\(/ + d\ir* + zbsexp{iilt) , 



£2m 



1 d^u 



ld\' 



2dz^ 

which, in turn, may be reduced to 



(18) 



(19) 



Q?-—(u-v) = 
oz 



'2dz^ 



1 -1 



— (m — v) — £1^5 . 
oz 



(20) 

The solution to Eq. ( |20] | is a superposition of a running 
wave Ui exp (/£2f — iqz) and spatially uniform oscillations 
Uoexp{iQ.t), 



dz 



Uo + Ui exp{-iqz), 



(21) 



where Oq i are the respective amplitudes. The running wave is 
the Bogoliubov elementary excitation (the quantum acoustic 
wave) with the wave number q{Q.) determined via the Bogoli- 
ubov spectrum 



(22) 



It follows from the boundary condition at the wall that Uq 

-Ui = -bjQ. and 



5f/ = ^ [sin {Q.t - qz) - sin£2f] . 



(23) 



5 



The perturbations of density follow from the ID continuity 
equation 

dt dz 



5n = ^ sin {£lt - qz) ■ 



(24) 



Then, the condition of weak perturbations 5n <C no becomes 
or in the dimensional form 



2mD? 



«1, 



(25) 



where q and Ci are the dimensional wave number and fre- 
quency of the induced excitations, coupled by the Bogoliubov 
spectrum (l22T i. In the low energy limit q"^ <C q'^, the Boguli- 
ubov excitations are sound waves with velocity ~ ^Jgho/m 
(unity in our dimensionless units). Then the restriction on the 
amplitude of the magnetic force reads 



(26) 



The criterion ( |26] | may be interpreted as the condition of 
incompressible condensates with the speed of sound much 
larger than any other parameter of velocity dimension in- 
volved in the problem. A more general condition dZSl ) takes 
into account the Bogoliubov spectrum for higher energies. 

We solved the Bogoliubov equations for a model infinite 
system. The case of a finite system, confined by two walls, is 
solved analogously, and the boundary condition ( [TtI i should 
be imposed for both walls. Real systems are always finite, 
and therefore the maximal wavelength is limited by the size 
of the system Ly. In the low-frequency limit condition ( |25l l 
is equivalent to the energetic condition (|4|i. In this limit, the 
Bogoliubov excitations are sound waves, with maximal wave- 
length of order Ly\ substituting — Inq^^ to Q, and using 
the sound spectrum, one obtains i 



IV. LINEAR ANALYSIS OF THE INTERFACE DYNAMICS 

We now consider the dynamics of the interface between two 
infinite BECs in the presence of an external force, which may 
be time-dependent, b{t). Similarly to Refs. jsi [33i [ssll we 
use a variational approach. We represent the wave functions 
in terms of density and velocity potential as y/,' = 
■^/Tfjexp so that the action for the system of two BECs is 



5 = 

















^) - 








1 ( 2 




2(«i 



50,- , 1 



dt 



-zb{t) {n2-n\ 



(27) 



For the variational calculation, we employ an ansatz that is 
able to describe bending and compression of the interface 
between the two condensates. We are interested in pertur- 
bation wavelengths much larger than the interface thickness, 
kZp = kj ^JTy 'SC. 1, so that the density deformations may be 
treated as a bending of the unperturbed profile as a whole, 

nj{z,x,t) =njo[z-5zj{x,t)], (28) 

where njo are the stationary density profiles for the BECs. The 
variational functions are taken to be of the form of Fourier 
modes with separation of space and time dependence. 



dzj{x,t) — l^j{t)cos{kx), 



(l)jiz,x,t) = ^j{t)fj{z)cos{kx). 



(29) 



(30) 



In classical hydrodynamics, perturbations of both components 
of an inert interface are related by the continuity condition 
^1 = 1^2- In the quantum gas-dynamics of two BECs this is 
not necessarily the case, since we may have interpenetration 
of the wave functions. 

The Lagrange equations with respect to the variables S,j are 



Cj / dzn'jfj = £,j / dz 



fj^^i'^jfj)+k'fjnj 



(31) 



These Lagrange equations are satisfied if the respective conti- 
nuity equations are 



(32) 



First, the time-dependent terms on each side have to be 
equated, which, up to constant factors, leads to 



0(0-^.(0- 



(33) 



This reflects the fact that and are collective canonically 
conjugate variables. For the functions fj, Ref.|[^ gave the 
approximate solution 



Mz) 



if 
if 



z<0, 
z>0. 



(34) 



The same solution holds in the presence of a small external 
force under the conditions specified in Sec. III. We take Eq. 
( |34| ) as our ansatz for fj. More elaborate expressions for f\ [z), 
and fi [z] may provide a smooth transition in the interface re- 
gion; however, the ansatz Eq. ( l34b yields reliable analytical 
results, as we will demonstrate below by comparison to the 
numerical solution to the problem. 

We now derive the Lagrange equation for No time 

derivatives in appear in the Lagrangian, so the Lagrange 
equation reads 



0: 



5S_ 
5C; 



(35) 



where we introduce //(^j to label the nonvanishing terms of 
5S/5^j; the numbering I in //(yj refer to the respective terms 
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in Eq. dZTl ). For convenience, we assume a finite size of the 
system, L^, L^, though the final resuh does not depend on the 
size restrictions. We have 

dzdxnQj{z + 5zj)fj^jCos^{kx) = 



dzn){z)fj{z) 



(36) 



The next nonvanishing term comes from the zero-point mo- 
tion energy. 



1_5 

2 5C 



dx 



dz 



(37) 



The first term within the parentheses yields a constant inde- 
pendent of C,, and the second term equals 



2 2 5^1 



dz 



dx 

2 



^b{t)j^^ j dz![z!n,2{z!)~5z2n,2{7!)] 



The integrals can be calculated as 

Lr. , . 5 ( 1 



'5 1 



2 ' '5Ci V 2 



^(Ot7-( + --yKOCi, (40) 



b{t)C2- (41) 



The next term of 55/5 does not contribute to the equations 
of motion, while the last (seventh) term Ij is important. Anal- 
ysis of the terms calculated so far shows that the variables 
are governed by two coupled oscillator equations, and the cou- 
pling happens because of /?. Since the equations for and ^2 
are identical, the normal modes may be in-phase, ^1 = ^2, 
and counter-phase, = —^2- First, let us calculate Ij for 
the case of ^1 = ^2, which we denote by /7(,„). In this case 
5zi = 5z2 = 5z, and we find 



n{m) 



- / <i^r(l + 7)nin2 



L, 5 
T50 



J dz (1+7) «oi (z + dz) «02 (z + Sz) = 0. (42) 



Next we consider the case of counter-phase oscillations with 

5zi = — 5z2 = 5z and ^1 = — ^2 = C' which leads to 



(38) 



The next non-vanishing term corresponds to the perturbation 
of the potential energy of the system. 



^5(;) ^J^. I ^ '■^^(^) ("2 - til) 

hi 
2 

L 5 f 

y(-l)'^(f)g^ J dzznoj{z+8zj). (39) 

It is more convenient to calculate these terms for /' = 1 and 
for j = 2 separately. Using the coordinate transformation z' — 
z + 5zj we modify the last expression as 



h(co) = y ^ / ( 1 + ^) «oi {z) «02 (z' - 25z) . 

To calculate 5/5^ in the last expression, it is easiest to as- 
sume 5z small and Taylor expand no2 {z! — 25z) to second or- 
der, which results in 

^7M = y4(l+7)C j dznoi-^. 



The last result may be transformed to a symmetric form simi- 
lar to 111: 



/7(.„) = -y4(l + 7)Cy^z — — . (43) 



L 5 f 

^5(1) = J dzznQi{z + 5zi) 



Taking into account all the terms obtained above we write 
equations of motion for ^j. The integrals are calculated us- 
ing the density profile ( fTTT ) 



dz' [z'noi {z') ~ 5zinoi {z')] 



-Lz+Szi 



(44) 



^5(2) 



J dzznQ2 (z + 5z2) 



dz I ^\/"o> 



V8, 



(45) 



7 



J 



dz— — = -V87/3- 

dz dz 



(46) 



We designate the interface perturbations of the in-phase oscil- 
lations by = ^2 = Cin ™d obtain the final equation 



5f2 



^-bit)k 



0. 



(47) ^ 



In the opposite case of the counter-phase perturbations = 
-C2 = we find 



^'a/|+^(1 + 7)V87-K0 



k = Q. (48) 



In the limiting case ofb = our theory reproduces the results 
of lH]. Within the hmit of a weak field, b <C b,n = J-J^, and 
long wavelength perturbations k/^/2y<^ 1, the terms in the 
parenthesis of Eq. ( |48l ) are of different orders of magnitude, 
with the term 4(1 + 7) \/87/3 dominating. For this reason, to 
the leading terms, the counter-phase mode is not affected by 
the external potential, and we do not consider it in the present 
work. In the following, we work only with the in-phase mode 
omitting the subscript "in" as C, [t) = Qn{t), and use Eq. (|47] | 
as the starting point for the subsequent analysis. 



V. DYNAMICS UNDER DIFFERENT TYPES OF 
TIME-DEPENDENT DRIVING FORCE 

In this section we demonstrate the diversity of hydrody- 
namic effects in a system of two separated BECs under the 
action of a time-dependent force. In particular, we are inter- 
ested in the RT instability developing due to a constant force, 
in the RM instability for a pulse force, in the parametric in- 
stability for an oscillating force, and in the interface dynamics 
under the action of a stochastic Gaussian force. 



A. The RT instability driven by a constant force 

Here we consider a constant force b (f ) =bQ= const, which 
drives the RT instability and capillary-gravitational waves in 
the unstable or stable configurations, respectively. Constant 
coefficients of Eq. (|47| i lead to harmonic oscillations of the 
interface ^ (f ) oc exp(/(Bf ) with the dispersion relation 



(0 



:k'^-kb0, 



which may be written in dimensional form as 



6? 



-ih n 



2rrP- 



^^2n(ax2-a)-k- 



2m 



(49) 



(50) 



where (b and k are the dimensional perturbation frequency and 
wave number. Note that Eq. (|49] | has the same form as the 
dispersion relation for capillary-gravitational waves in clas- 
sical hydrodynamics Ii34i] . When the external magnetic field 
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FIG. 3: Lowest Bogoliubov modes (the RT dispersion relation) for 
different values of the driving force. The analytical result is obtained 
from Eq. l |49l l, and presented by lines. Dashed lines show real parts 
of the Bogoliubov eigenvalues (waves), and solid lines correspond to 
the imaginary parts (instability). The markers show respective nu- 
merical solution to the linearized GP equations; filled markers cor- 
respond to the instability regime and empty markers correspond to 
waves. 



is zero, we find capillary waves with frequency related to the 
wave number as 



al^k\ 



(51) 



In the case of positive gradient of the magnetic field the dis- 
persion relation Eq. ( |49] l yields < for perturbations with 
sufficiently long wavelengths, which corresponds to the RT 
instability ^ (f) oc exp(aoO with the growth rate depending on 
the wavelength as 



an 



-bQk{\-k^/k]) 



and with the cut-off wave number 



(52) 



(53) 



According to Eq. ( |49l ), the maximal growth rate corresponds 
to k,„c,x = kc/ \/3- Figure 3 compares the analytical results of 
Eq. (l49b . ( fSOb to the numerical solution to the linearized GP 
equations. The numerical solution is obtained by diagonal- 
ization of the linearized GP equations, approximated by finite 
difference method on a grid of 3335 points on the interval 
[—240; 240] in the dimensionless coordinates. Figure 3 shows 
good agreement of the analytical theory with the numerical 
results. The instability growth rate Eq. ( |52] | agrees also with 
the formulas employed in lfT2ll . 

Though the present work considers only the linear stage of 
the RT instability, for illustrative purposes, it is worth present- 
ing the results obtained for the nonlinear stage. For exam- 
ple. Fig. 4 shows the development of the RT instability at the 
nonlinear stage in a confined system of two BECs obtained 
numerically using the methods of jlSll . In the numerical so- 
lution, we begin with a system of 3.2- 10^ atoms of **^Rb, 
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FIG. 4: Development of the RT instability in a trapped system of 
two-component BEC presented for one of the components. 
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equally split between the two states, with |F,mf) = 1) 
and 1 1,1). We use the "pancake" trap geometry with (B, ~ 
m^ = 2n- lOOHz and — 1% ■ 5kHz. We solve the coupled 
GP equations numerically to find the ground state of the sys- 
tem. The two components occupy the regions with z > and 
z < 0, respectively. A perturbation of size 9/im is used, and 
the field of view for all images is 55/xm square. The mag- 
netic field of magnitude B' — 1.78G/m is added to the sys- 
tem, directed so that the two components are pushed toward 
each other The system is shown at times 15.9, 19.08, 22.26, 
25.44 (ms) after the magnetic gradient is added. The nonlin- 
ear stage may be described qualitatively as development of the 
initial sinusoidal perturbations into a mushroom structure with 
subsequent generation of quantum vortices. More detailed in- 
vestigation of the nonlinear stage of the RT instability will 
be presented elsewhere. Similar results may be found also in 



be nrt 

ill. 



B. High-frequency stabilization of the RT instability due to an 
oscillating force and parametric instability 

In this subsection we consider an oscillating force h [t) — 
bo + bsCos{Q.t), where bo and bs are some constant ampli- 
tudes and Q. is the frequency of the force. We show that an 
oscillating force leads to two interesting effects: 1) High- 
frequency oscillations may stabilize the RT instability pro- 
duced by the constant background force bo; 2) The oscillat- 
ing force may trigger a parametric instability even in the ab- 
sence of any constant acceleration, e.g. for bo — 0. We start 
with the first effect, which is similar to the mechanical phe- 
nomenon known as the Kapitsa pendulum fs^. We assume 
that the frequency Q. of the external force is much larger than 
the RT instability growth rate Oo = y^bok{l — k^/k^), that is, 
£1^ (X(). Then we separate the slow averaged terms (labeled 
"a") and fast oscillating terms (labeled "f ) in the solution to 
Eq. (|47]i as C = C« (0 + C/ (0 cos (Of), where Ca (f ) and C/ (f ) 
change slowly on the time scale of high-frequency oscillations 
2n/Q.. Following the calculation method of |30], we take into 
account the relation ^ w Q, — Q?l^f cos {Q.t) and modify Eq. 
(gill to 

C„ - ii^C/cos (at) = bok (i - ^"j [c„ + C/cos {at)] 



FIG. 5: Stabilization of the RT instabiUty by a high-frequency force 
for the ampUtude ratios bs/bo = 0;7; 15;20, for bo = 2.5 ■ 10"^, k = 



-bsk 



Ca COS {at) + ^ C/ + ^ C/COS (2i2f ) 



Separating the terms oscillating with different frequencies we 
find the relation between the amplitudes of the slow and fast 
terms 

„ bsk „ bgk „ 

" ~a2 + bok{l-k2/k^J^" " £22 ^" 



and reduce Eq. (l47b to 



C„ = bok[l 



k^ bik 



(54) 



Equation ( |54] | describes stabilization of the RT instability by 
high-frequency oscillations for 



k^ b]k 
k^.~2b^ 



<0. 



Thus, stabilization of the RT perturbations of wave number 
k by the high-frequency term is expected for the oscillation 
amplitude 



(55) 




which may be written in dimensional units as 



B' 



> 



PLbB'Jc 



(56) 



The stabilization of the RT instability by high-frequency os- 
cillations is illustrated in Figure 5, which presents the results 
of direct numerical simulation of Eq. Wf\ for the amplitude 
ratios bs/bo = 0;7; 15; 20, for the stationary component of the 
magnetic field bo = 2.5 • 10^^, the wave number correspond- 
ing to the maximal growth rate k = k„„,x and frequency of the 
oscillating magnetic field component £2 = 20ao. In agreement 
with the theoretical predictions, unstable exponential growth 
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FTG. 6: Stability limits (white line) in the parameter space of the 
scaled wave number k/kc and the oscillation amplitude bs/bo, where 
bQ = 2.5 • 10^^, £1 = 20max(o!()). Shading indicates value of the 
instability growth rate in the unstable region, Eq. 
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FIG. 7: The RT instability growth rate, Eq. l l54b . in presence of 
external high-frequency field with fl = 20max(ao) vs perturbation 
wave number for different values of b\/b() — 0;7; 15;20 and b{) = 
2.5 -10-^ 



of perturbations (in average) for hsjh^ — 0;7 is replaced by 
waves at the interface for sufficiently large amplitudes of the 
oscillating part of the magnetic field, bs/bo — 15; 20. Figure 
6 shows stability limits in the parameter space of the scaled 
wave number k/kc and the oscillation amplitude b^/bo; the 
shading indicates the absolute value of the instability growth 
rate in the unstable region. Figure 7 shows the RT instabil- 
ity growth rate in the presence of the external high-frequency 
magnetic field plotted versus the perturbation wave number 
for different values of the external field. 

The second effect related to the oscillating external field 
is the parametric instability, which is known from mechanics 
ll36h and which has been encountered in gas dynamics, e.g., 
of combustion systems |30- 32]. The constant part of the "ac- 
celeration" is not required to obtain the parametric instability 
and it may be omitted in the calculations by taking bo = 0. 



Similar to ll30l 13611 we look for a solution to Eq. ( |47] | in the 
form = [^i cos (Hf /2) + ^2 sin (iif /2)] exp (at). Mark that 
the growth of the parametric instability is accompanied by os- 
cillations of the perturbation with frequency fi/2, i.e. half of 
the driving force frequency. Substituting ^ into Eq. (l47b we 
obtain 



n2 
4 



£,1 cos 



: bsk COS 



+OI0Q. 

at 



^1 cos 



h E,2 sin 
-(^2 cos 
& sin 



T 



Separating the terms with cos (i2r/2) and sin(i2r/2) and omit- 
ting higher-frequency terms, we reduce the equation to the lin- 
ear algebraic system for the amplitudes and ^2 



Oo 



4 



Solving the system we obtain the dispersion relation 



2 a 



bik 



a^Q} = 0. (57) 



Substituting a = into Eq. dSTl i we find the limiting oscilla- 
tion amplitude b^ required to induce the parametric instability 



2 £22 
Or 



In dimensional variables the last relation reads as 

^>k-'\Aal-Q?\, 

m ' ' 

where the dimensional frequency of capillary waves is 



(58) 



(59) 



£2, [k] 



2-1/2 







2m2 



Minimal (infinitely small) critical amplitude of the magnetic 
field corresponds to the perturbations with the wave number 
determined by the equation £2^.(A:) — £2/2. Therefore, when 
applying an external magnetic force b (f) = bscos (£2f) to the 
system of two BECs, we should expect resonance pumping of 
the capillary waves with frequency £2/2 and the wave number 
A: = £22/3(27)1/6. If capillary waves of another non-resonant 
wave number dominate initially, then pumping of the waves 
requires finite (not necessarily small) amplitude of the driving 
force. The instability domain is shown in Fig. 8 in the param- 
eter space of the force amplitude and the wave number with 
the shading indicating the perturbation growth rate. Figure 9 
presents the dispersion relation for the parametric instability 
(the growth rate versus the wave number) for different ampli- 
tude values of the oscillating magnetic field. 
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which yields 



FIG. 8: Stability limits (white line) in the parameter space of the 
wave number k and amplitude bs of the external high-frequency force 
with fl = 3.6- 10^-'. Shading indicates value of the instability growth 
rate in the unstable region, obtained from Eq. ([57}. 
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FIG. 9: Growth rate of the parametric instability, Eq. J57b . as 
function of k for different amplitudes of the external oscillating 
magnetic force = 0.25 • 10""*; 1.75 ■ 10""*; 3.5 ■ lO"'*; 5 • lO"'* with 
n = 3.6- 10"^ 



C. RM instability produced by a magnetic pulse 

Within the context of quantum gases, the RM instability 
produced by a magnetic pulse b{t) — b^iS (f) has been dis- 
cussed recently in Ref. |15], which indicated considerable 
differences in the instability development in BECs in compar- 
ison with the traditional configuration of classical gases sep- 
arated by an inert interface. Here we analyze the linear inter- 
face dynamics by a more general (matrix) method of solving 
Eq. ( |49] l applicable to a wide class of time-dependent forces. 
Particularly, the matrix method will be extended to the case 
of stochastic forces in the next subsection. We introduce the 
vector function = ( ^ ^ ) and rewrite Eq. ( |49] | as 







(60) 



where the matrix may be split into a constant and time- 
dependent parts. First, Eq. ( l60t is solved for the case b — 0, 



1 

-a? 



Cvit)= exp 
cos (fie?) Q.^^sm{Q.ct 



-D.cSm{D.ct) cos(i2tf) 
For convenience, we introduce the designations 



C..(o) 
Cv(o). 





1 



(61) 



(62) 



Gow^f ^^";A^ ^^'''<^fA. (63) 

^'^ ' \^-i2(.sin(ntf) cos(ntf) J ' ' 

We look for a general solution to Eq. ( l60l l in the form (0 = 
Go (t) V (t) and reduce Eq. ( l60b to 



v{t) = kb{t)S(t)v(t), 



(64) 



where 



V(f) = Go^(f)VoGo(f) 



(2^^ sin (2i2,.f ) -D.-^'sax?- (il^f ) 
cos2 (act) (2£2cr' sin (2£2,.f 

The solution to Eq.(l64b is the time-ordered exponential F (t), 
e.g. see [,37,] . 



exp 



I 



dTF(t) 



Y.JdhJdti--- |^/r„F(fi)F(f2)---F(f„). (65) 







We solve Eq. ( |64l l, go back to the original variable and obtain 
the general solution to Eq. ( |60l ) 



Cv (f) = Go (f)exp J jdskb{s)y{s) \ Cv(0). (66) 



On the basis of the general solution, we may recover the re- 
sults of Ref. ifTsIl on the RM instability triggered by a pulse 
b{t) = b^i5{t). In this case all terms of time-ordered exponen- 
tial commute with each other and we obtain 



exp< J dskbcid{s)\{s) 



- exp 



r 








[kbd 




[kb, 
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FIG. 10: Development of the RM instability in a trapped system of 
two-component BEC presented for one of the components. 



which yields the final solution 



at) - J a (0) + ^ [kb,C (0) + C (0)J sin {a,t + <p,) 



with the new oscillation phase determined by 

a.C(o) 



tan <P(i 



ao)+kb,t;{o) 



(67) 



(68) 



Thus, the general solution Eq. ( |66] | reproduces the analytical 
results of Ref. ifisll given by Eq. (l67l i. Particularly, Eq. ( l67l i 
describes the jump in the amplitude and phase of a standing 
capillary wave because of the pulse for a non-zero bending of 
the interface, (0) ^ 0. In the opposite case of (0) = and 
1^ (0) 7^ the pulse does not influence the interface dynamics. 
Eq. (l67T i shows that the transfer of magnetic energy of the 
pulse to superfluid oscillations is the most efficient when the 
deviation of the interface just before the pulse is maximal. 

The linear analysis works rigorously as long as the pertur- 
bation amplitude after the pulse remains sufficiently small. 
Once it becomes large enough, nonlinear effects dominate, 
and the interface evolves into a mushroom pattern with de- 
tachment of droplets as discussed in |15]. Figure 10 illus- 
trates the development of the RM instability in BECs at the 
nonlinear stage obtained using the numerical methods similar 
to ifTsll . We use the same system of BECs as in Fig. 4. A mag- 
netic pulse of duration T = 0.02ms and size y/nnsB'T/lm = 
1.05 • lO^^m/s is applied to the initial capillary wave seen in 
Fig. 10 (a). The following graphs correspond to the time in- 
stants 2.39; 4.77; and 7.16 ms after the pulse (figures b, c, d, 
respectively). Development of the RM instability at the initial 
stages (e.g. Fig. 10 b) resembles the RT instability shown in 
Fig. 4. However, the difference between two instabilities is 
noticeable at later stages shown in Fig. 10 (c, d). In particular, 
we observe the detachment of droplets in the RM instability 
in agreement with the model presented in Ref. lITsIl . 



D. Gaussian noise 

Finally, we consider the interface dynamics under the action 
of a Gaussian random force br (f ) with mean zero {b,- (f )) — 0, 
where the brackets imply averaging over a statistical ensemble 



of realizations. The force is switched on at f = 0. Because of 
the Gaussian property, the force is characterized only by vari- 
ance Q and by the coiTelation time T^. The case of zero is 
known as the Markovian case, and in this case averaging may 
be performed exactly. On the contrary, the case of non-zero 
T( is not solvable in general, but an analytical approximate so- 
lution may be obtained in the limit of small correlation time 

In order to obtain an averaged form of Eq. ( |47] ) we em- 
ploy the general solution ( |66] |. Since Eq. ( |66] | is a time- 
ordered exponential with non-commutative kernel at different 
moments of time, we use the cumulant method for averaging 
iIstIi . Within this method, if G(t) is a stochastic Gaussian 
operator with mean zero (G(r)) — 0, then the time-ordered 
exponential of G(t) is averaged by the formula 



exp 



j dTG{T) 


^ =exp 


.0 





dx / ds{G{T)G{s)) 



(69) 

In the Markovian case this formula is exact, but in the non- 
Markovian case this is an approximation valid when is 
small compared to the intrinsic time scale of the system, 
Ttilc <^ 1, and the noise force is sufficiently weak to be treated 
as a perturbation |37]. The kernel of the right-hand side of 
Eq.(|69ll is the cumulant series terminated at the first non- 
vanishing term (the first-order term vanishes because (G(f)) = 
0). We note that all cumulants in Eq. (6% vanish identically 
for Tc —0 because of the Gaussian property, but they are non- 
zero for finite due to non-commutativity of the kernel in Eq. 



We start with the Markovian case when the autocoiTelation 
function is given by the delta-function 



br{t)br{s))=Q5{t-s). 



(70) 



Using Eqs. ( |66l l, ( |69] l, we find that the averaged solution to 
Eq. ( l47l i reads 

(C(0> = C(0)cos + ^ sin (n,f) , (71) 

i.e. the Markovian stochastic force does not influence aver- 
aged dynamics of the interface. 

Now we consider a non-Markovian exponentially corre- 
lated noise with autocorrelation function given by 

{br{t)br{s)) = :^exp (- \t - s\ /T,) , (72) 

which reduces to the Markovian case in the limit of zero cor- 
relation time T( 0. We average Eq. (l47l l using the Novikov 
theorem L38il : 



d^ 
dt^ 



m)+^lm)+k j ds{br{t)br{s)) 



5b, .(s) 



(73) 
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We calculate the functional derivative in Eq.(l73Tl using the 
general solution Eq.( 



Go(f)exp<^ / dskb{s)\{s) } Cy(0) 



5b, (s) dbr (s) ■ 
- dit-s)9is)Goit)Go\s)k\QGo{s) 

X exp J / dskb{s)Y{s) \ Cv(0) +A[/7,] (f,i)Cv(0) 



. (f - (.)Go {t - ^) WoGo \t-s) Cv(f ) 
-A[br] (r,^)Cv(0), 



(74) 



where 9{t) is the Heaviside step function, and A[bi-] (t.s) is 
a functional arising from the non-commutativity of the kernel 
of the exponential in Eq.(|74li. The functional A [br] {t,s) may 
be dropped in the present limit of weak noise. We obtain from 
Eq. (Elll 



5br {s) 



d{t~s)d{s)k{[y{t~s)],,m)- 



[v(f-^)]i2^(C(0> 



(75) 



Substituting the operators Go and Vq to Eq. ( [75] ). and using 
Eq. ( |73] l we find 



^(C(f)) + n2(C(0)- 



/ 



ds—{br{t)b,{s))^m2£l,{t~s){l;{t))- 



2 

J ds^{brit)b,-is))sm^n,{t-s)j^{at))=0. 


(76) 

Calculating the integrals on time scales much larger than the 
correlation time of the noise f Tc, we amve to an equation 
describing the averaged interface dynamics 



dt^ 
where 



^ (C(f )) - 2^1 (C(f )) + (il? + v2) m) = 0, (77) 



1 k^QT^ 

/f - — 



2 1 + (2i2,T,)2 
is the average growth rate and 

1 k^QTc 



2 1 + (2n,T,) 



(78) 



(79) 



modifies the frequency due to noise. The solution to Eq. (fTTT i 
grows exponentially in time as 



^ o<: exp 



(80) 




FIG. 11: (a) Single realization of the external stochastic Gaussian 
Markovian force with = 5 • 10"^, = 50 and flcTc = 1.5 • 10"^. 
(b) The respective numerical solution to Eq. l |47t . The dashed line 
shows amplitude of the averaged analytical solution, Eq. J80t . 
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FIG. 12: (a) Single realization of the external stochastic Gaussian 
Markovian force with fi = 2 ■ 10"^ and = 1200 and flcfc = 0.38. 
(b) The respective numerical solution to Eq. l l47t (two realizations). 
The dashed line shows amplitude of the averaged analytical solution, 
Eq. liOll. 



Thus, a stochastic non-Markovian force leads to unstable 
growth of the average perturbation amplitude with modified 
frequency of the capillary oscillations. The analytical theory 
for interface dynamics under the action of a stochastic force 
is compared in Figs. 11, 12 to the numerical solution to Eq. 
(|47] l calculated for k = 0.0154. We note that real stochas- 
tic processes always have finite frequency band width, though 
it may be very large. Therefore, any real stochastic process 
is characterized by some non-zero correlation time leading 
to exponential growth of amplitude of the oscillator Eq.d47l), 
which means that the Markovian case may be recovered only 
as a limit of ilcT; 0. Figure 11 shows the numerical so- 
lution to Eq. (4T\ for the stochastic force with Q — 5 ■ 10^^ 
and ClcTc = 1.5 • 10^^ <C 1, which may be considered as a 
good approximation for the Markovian case. The dashed line 
represents the envelope function of the averaged amplitude of 
the interface oscillations predicted by the theory, Eqs. dTSl ). 
dSOl ). In the Markovian case, the averaged amplitude does not 
change, which is supported by the numerical results shown in 
Fig. 11. Figure 12 corresponds to a non-Markovian case with 
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Q = 2- IQ-'^, Zc = 1200 and H^T,. = 0.38. The numerical solu- 
tion to Eq. ( |47] l is shown by the solid lines for two realizations 
of the noise, while the dashed line represents envelope of the 
averaged amplitude increasing exponentially in time accord- 
ing to the analytical theory Eqs. ( iTSl l. (ISOl l. In agreement with 
the theoretical predictions, the numerical solution also demon- 
strates exponential increase of the perturbation amplitude. 

VI. SUMMARY 

We have studied the dynamics of an interface in a two- 
component BEC driven by a spatially uniform time-dependent 
force. Experimentally, such a force may be produced by a 
magnetic field gradient acting on the system of two interact- 
ing BECs with different spins |fl2|. By applying the varia- 
tional principle to the GP Lagrangian, we have derived the dis- 
persion relation (or the respective ordinary differential equa- 
tions) for linear waves and instabilities at the interface. We 
have considered different time-dependent forces leading to a 
diverse collection of dynamical effects: (i) A constant force 
pushing the components towards each other generates the RT 
instability for a certain domain of perturbation wave numbers. 
Because of the RT instability, small perturbations grow expo- 
nentially in time without oscillations, (ii) A sinusoidal force 
with frequency £2 leads to the parametric instability at the in- 



terface. The critical strength of the force required to drive 
the instability depends on the perturbation wave number. This 
strength is infinitely small for the intrinsic capillary waves at 
the interface with frequency equal to £2/2. Because of the 
parametric instability, small perturbations not only grow expo- 
nentially in time, but also oscillate with frequency fi/2. (iii) 
A pulse force leads to the quantum counterpart of the RM in- 
stability. Within the linear approximation, the force produces 
a discontinuous jump in the amplitude and phase of the cap- 
illary waves at the interface, (iv) A non-Markovian stochas- 
tic external field (with non-zero correlation time) drives the 
instability at the interface accompanied by oscillations. The 
growth rate of the instability is determined by the variance of 
the driving force and by the correlation time. The oscillation 
frequency is shifted in comparison to the intrinsic frequency 
of the capillary waves. A Markovian force on average does 
not lead to any effect. 
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